In [1]:
%matplotlib inline
pylab.rcParams['figure.figsize'] = (10.0, 8.0)
from numpy import *
from scipy.optimize import root
In [2]:
def population(p0, m, n, t):
return m*p0 * exp(m*t) / (n*p0*exp(m*t) + m - n*p0)
def ddm_population(p0, m, n, t):
"""Dérivée partielle par rapport à m
"""
return p0*exp(m*t)*(n*p0*(exp(m*t)-1)+t*m*(m-n*p0))/(n*p0*exp(m*t)+m-n*p0)**2
def ddn_population(p0, m, n, t):
"""Dérivée partielle par rapport à n
"""
return m*p0*p0*exp(m*t)*(exp(m*t)-1)/(n*p0*exp(m*t) + m - n*p0)**2
In [3]:
def solveur(p0, t1, p1, t2, p2, facteur=1, method="hybr"):
def residu(x):
"""Fonction donc on cherche le zéro
"""
m, n = x[0], x[1]
return array([
population(p0, m, n, t1) - p1,
population(p0, m, n, t2) - p2
])
def jacobien(x):
m, n = x[0], x[1]
return array([
[ddm_population(p0, m, n, t1), ddm_population(p0, m, n, t2)],
[ddn_population(p0, m, n, t1), ddn_population(p0, m, n, t2)]
])
def solveur_parametre(x0):
resultat = root(residu, x0, jac=jacobien, method=method)
return resultat
return solveur_parametre
In [4]:
pop_0 = 29981336 # Population initiale
pop_7 = 31359340 # Population après 7 ans
pop_13 = 32402940 # Population après 13 ans
Numériquement, le problème n'est pas très bien posé. On utilise donc le changement de paramètres suivant $$ p_{p_0, m, n} = \frac{1}{k}p_{k\, p_0, m, \frac{n}{k}}$$
In [7]:
facteur = 1e-7
x0 = array([1e-5, 1e-12/facteur]) # Point initial
s = solveur(pop_0*facteur, 7, pop_7*facteur, 13, pop_13*facteur, method="hybr")
resultat = s(x0)
m_sol = resultat.x[0]
n_sol = resultat.x[1]*facteur
print "m : {0}\nn : {1}\npopulation limite : {2}".format(m_sol, n_sol, m_sol/n_sol)
In [8]:
population_parametree = vectorize(lambda t:population(pop_0, m_sol, n_sol, t))
t = range(300)
plot(t, population_parametree(t))
plot([0, 7, 13], [pop_0, pop_7, pop_13], "o")
xlabel(u"Années")
ylabel("Population")
axhline(m_sol/n_sol)
Out[8]:
In [ ]: